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The past decade has seen a significant growth in research targeted at space 
based observatories for imaging exo-solar planets. The challenge is in design- 
ing an imaging system for high-contrast. Even with a perfect coronagraph 
that modifies the point spread function to achieve high-contrast, wavefront 
sensing and control is needed to correct the errors in the optics and generate a 
''dark hole" . The high-contrast imaging laboratory at Princeton University is 
equipped with two Boston Micromachines Kilo-DMs. We review here an algo- 
rithm designed to achieve high-contrast on both sides of the image plane while 
minimizing the stroke necessary from each deformable mirror (DM). This al- 
gorithm uses the first DM to correct for amplitude aberrations and the second 
DM to create a fiat wavefront in the pupil plane. We then show the first results 
obtained at Princeton with this correction algorithm, and we demonstrate a 
symmetric dark hole in monochromatic light. 

(c) 2011 Optical Society of America 
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1. Introduction 

With the need to image faint exo-planets close to their parent star, wavefront control for 
space-based high contrast imaging has received tremendous attention over the past few years. 
The idea of completely canceling starlight using coherent subtraction via a Deformable 
Mirror (DM) was first introduced by Malbet et al. [1], where they proved, using non linear 
numerical tools, the theoretical feasibility of a very high contrast Dark Hole in the image 
plane of a telescope. Experimentally the first idea tested was based on a linearization of 
this previous result: cycle through a set of arbitrary DM configurations and choose the one 
yielding the best contrast in an algorithm called "Speckle Nulling" (e.g.. Brown & Burrows 
[2], Trauger [3]). Borde & Traub [4] proposed a refinement of this solution that yielded 
much faster convergence rates: they separated the estimation and the correction stages 
and reduced the second one to a simple matrix inversion based on an energy minimization 
criteria. However this method presented numerical caveats when associated to a coronagraph 
since it required the inversion of an ill-conditioned matrix. Give'on et al. [6] proposed a 
solution to regularize the problem based on Electrical Field Conjugation. In this paper we 
introduce an alternative correction method that fully solves the non-linear problem while 
using the smallest DM deformation possible. We not only present the theory underlying 
these algorithms but also provide an experimental validation of the Stroke Minimization 
method using one and two DMs. Previous experimental results of high contrast Dark Hole 
formation only achieved high contrast in one-half of the image plane. Because of the unique 
dual DM feature of the Princeton University High Contrast Laboratory, this paper is the 
first experimental report of symmetric high contrast point spread functions. 
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One of the characteristics of high contrast wavefront control in space is focal plane wave- 
front sensing. Using the science camera to also estimate the wavefront avoids potential 
non-common path errors from conventional pupil estimation. The control algorithms dis- 
cussed in this paper, therefore, all assume perfect estimates of the field in the region of the 
image plane where we seek high contrast. The experiment uses an estimation algorithm that 
retrieves the wavefront based on random DM diversity (Give'on et al. [5] and Give'on et al. 
[6]). That is, probes are set on the DM in order to cause interference with the aberrated 
field. A careful analysis of the resulting speckles lead to the retrieval of the complex field. 

For all the analysis and simulation of the various control approaches in this paper we 
assume perfect estimates of the focal plane field. 

2. The Optimal Dark Hole Problem 

The goal of the correction problem is to cancel the starlight due to optical imperfection so 
that the contrast in the image plane is the one for which the coronagraph was designed. It 
was first shown in Malbet et al. [1] that this could be formulated as a nonlinear optimization 
problem. While the problem itself is straightforward to state, many solution approaches are 
possible. Here, we follow the notation of Give'on et al. [6]. We model the coronagraph as 
a general linear transformation, (7, between the electric field at the deformable mirror, Eq^ 
and the final electric field at the image plane of the camera, 

Ef = C{Eo}. (1) 

When the DM is at the pupil, the linear transformation C is a Fourier transform. When the 
DM is not at the pupil, the operator includes Fresnel propagations from the DM to the pupil. 
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The linear mapping could represent any number of possible coronagraphs (shaped pupils, 
Lyot type coronagraphs, pupil mapping, etc.). In most cases, it will simply be a series of 
Fourier transforms and convolutions involving various mask types. For the simulations in 
this paper, we assume a shaped pupil coronagraph because of its simplicity. There, C{Eo} 
simply represents the Fourier transform of the field at the shaped pupil. For this analysis, 
the DM is assumed to be in a plane conjugate with the shaped pupil, while the experimental 
results presented in later sections, take into account propagation between the DMs and the 
plane of the shaped pupil. 

If we include amplitude and phase aberrations and a single DM to correct, the electric 
field at the DM can be written as (Give'on et al. [6]), 

where ax{x^y) and I3\{x^y) are respectively the amplitude and phase aberrations across 
the pupil (and may differ with wavelength due to propagation effects), A{x^y) is the pupil 
apodization, and ijj{x^y) is the DM height in units of wavefront. The problem is to find a 
DM surface to sufficiently cancel the aberrations and restore the contrast in some region of 
the image plane, often called the "dark hole" . 

This is an infinite-dimensional problem and cannot be solved exactly. One approach 
to making it tractable is to approximate the DM surface height by a finite sum of basis 
functions, 

N 

V^(x, 1/) = Ao X! CLkfk{x, y). (3) 
k=i 

where Aq is the central wavelength in the band considered. The problem is now to choose the 
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coefficients to achieve the desired performance. The basis functions can be any convenient 
set, such as Legendre polynomials, Zernike's, Fourier modes, Chebyshev polynomials, or 
others. Depending upon the choice of basis function, the coefficients, a^^ will be some 
function of the voltages applied to the DM. Selecting a particular expansion is a tradeoff 
among various objectives. One objective might be to minimize the number of terms for an 
adequate fit. Another might be to simplify the relationship between the coefficients and 
the voltages. Most often, the fk are chosen to be the so-called influence functions. These 
are the shapes the DM takes when a voltage is applied to only a single actuator. The 
assumption is then made that superposition holds (that is, that an arbitrary shape can be 
found by summing infiuence functions). For this choice, the coefficients have a one-to-one 
correspondence to the voltage applied to the k^^ actuator. The experimental work reported 
in this paper uses such an influence function basis. 

To further simplify our notation, we write the coefficients to be found as a single column 
matrix, X — [ai, a2, . . . , aAr-i, tiAr]^ and the basis functions in the row matrix F{x^y) — 
y), /2(x, y), . . . , /Ar-i(x, y), /Ar(^, 1/)], which lets us write the DM surface in matrix 

form. 



Finally, for convenience we streamline our notation by replacing the complex exponentials 
in Eq. (2) by a function /i. 



where 7a = + is the complex, wavelength dependent abberation and h represents the 
nonlinear dependence on the aberration and DM settings. 



(4) 



= A(x,j/)/i(7a(x,j/),X) 



(5) 
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The simplest correction algorithm would be to choose the coefficients such that the DM 
surface cancels, or conjugates^ the aberrations in Eq. (2). This is the typical approach used 
in conventional, ground based adaptive optics (AO). We notice immediately two limitations. 
First, with a single DM it is only possible to correct the field at a single wavelength. Second, 
since the DM is a phase correcting device, it nominally only corrects for the phase aberration. 
However, complete phase conjugation across the pupil is more than is needed, as we seek a 
dark hole in only a limited region of the image plane. In fact, we show in the appendix that 
by correcting only the part of y) whose image plane electrical field is centro-symmetric 
with respect to the optics axis , it is possible to correct for both amplitude and phase in a 
dark hole on only one half of the image plane with a single DM. For symmetric dark holes, 
two DMs are necessary. In Pueyo & Kasdin [9], we show how, in principle, two DMs can 
be used to correct for both amplitude and phase across the image plane and at multiple 
wavelengths. We discuss two DM corrections in Section 3. 

Deterministic phase conjugation also requires knowledge of /3a (x, y). In classical AO, this 
is done via a wavefront sensor in a diverted beam. For extreme high contrast in space, 
however, non-common path errors result in an uncorrectable component that degrades the 
desired contrast. Our goal is to develop algorithms for correction using only measurements 
of the image plane field, avoiding the need to backpropagate to estimates of the pupil 
aberrations. In the remainder of this section we review past approaches and present a 
new approach, stroke minimization. We also continue to restrict ourselves to algorithms 
incorporating one DM at a single wavelength. In section 3 we show how to generalize this 
one DM solution to two sequential DMs and present the first experimental validation of two 
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DM dark hole generation. 
A. Energy Minimization 

One approach to achieving a dark hole is to minimize the energy in the desired region of 
the image plane as first proposed by Malbet et al. [1]. While this first form of the problem 
resulted in a complicated non-linear optimization, Borde & Traub [4] presented a simplified, 
linear form which we review in this section. 

The total energy in some region S of the image plane is written as the integral of the 
intensity over that region, 

Ss = {Ef,Ef)^ = J J^EfE}d^dri, (6) 

where * represents the complex conjugate and (^, 77) are image plane coordinates and Ef — 
C{Eo}. 

In Energy Minimization we seek to minimize £s over X. There are a number of ways 
one might do this. For example, if we let g{X) — (Ef^Ef)^^ then a simple steepest descent 
algorithm gives, 

X„+i = x„ - pVxgiXr,) (7) 

for some constant z/, where Vxg{Xn) is the gradient of g{X) evaluated at X^. While simple 
to formulate, steepest descent is known to take many step to converge down a long narrow 
valley and tends to overshoot farther than the minimum (Press et al. [7]). Instead, we might 
choose a Newton algorithm based on a Taylor series expansion of g{X)^ 

Xn+l =Xn- H-^Vxg{Xn) (8) 
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where Hn is the hessian matrix of g{X) at X = X^. However, it is a common problem with 
Newton's method that when far from the solution, the Hessian can become negative definite 
and the step no longer descends. The usual approach is to use Gauss' correction, resulting 
in the Gauss-Newton iteration. 

The correction step to the coefficients is given by the pseudo-inverse of the Jacobian of the 
total energy in the dark hole. 

Using the expression for the electric field in Eqs. 1 and 2, the Jacobian can be written 

as: 



'c{Eo{X^)}, t'^C{Eo{X^)F}^ 



(10) 



where 3? stands for the real part of a complex number. Borde & Traub [4]'s insight was 
realizing that computing the gradient in the Gauss approximation to the Hessian around 
X = 0, and neglecting cross terms between ^x{x^y) and '0, would greatly simplify the 
problem at a small cost in convergence rate. Following their presentation, we write a first 
order expansion of h{jx{x^y)^X) as 

/i(7a(x, y),SXn) ^ 1 + 7a y) + y)5Xn (11) 

where we have expanded about X = and ignored the cross term in 
i^^F(x^y)5Xn 7A(x,y). Ignoring this cross term, leads to an errors of about |7a| per- 
cent in the actual gradient. This error in theory slows down the convergence rate. However, 
as discussed later in this section, for this application such an error in the gradient is not the 



convergence rate limiting factor. In any case, when the starting wavefront error is small, 
starting contrast below 10~^, the error on the gradient does not impact the minima towards 
which the algorithm converges. As a consequence we can write. 



'E,,^'-^C{AF}) 



(12) 



Vx^(X„)^Vx^(X„) ^ -2(2^)2 * 3? [{C{AF}, C{AF})s] • (13) 

This approximation significantly simplifies the algorithm as V xg{Xn)^V xg{Xn) can now 
be pre-computed and V xg{Xn) only depends on the projection of the image plane field at 
iteration n on the modal matrix C{AF{x^y)}. 

Implementing energy minimization thus requires only an estimate of the electric field 
at the image plane, Ef. Borde & Traub [4] suggested finding this field via image plane 
measurements only and a DM diversity estimation scheme, avoiding any non-common path 
errors associated with pupil sensors. They devised a reconstructor for Ej that used several 
DM settings to disentangle the ambiguity of intensity measurements at the science camera. 

Unfortunately, a common problem with the Gauss-Newton method is that 
V xg{XnYV xg{Xn) can drop rank. In fact, this is a signficant problem for dark hole gener- 
ation as the image plane area being minimized is small enough that is underdetermined. 
One common solution is the Levenberg algorithm, where a constant "damping parameter" 
is added to the search direction, 

x.+i = X. - {Vxg{Xr,fVxg{Xn) + ^l i)-'Vxg{Xnf- (14) 

where /i is a regularization parameter and / the identity matrix. This approach provides 
a balance between Gauss-Newton and steepest descent, allowing faster convergence when 
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far away and slower convergence (avoiding overshoot) as the minimum is approached. It 
also guarantees invertibility of the weighting matrix. Marquardt suggested a modification, 
replacing the identity matrix with the square of the Gauss approximation to the Hessian, 

This is known as the Levenberg-Marquardt algorithm. Marquardt 's damping correction 
is very similar to Tikhonov regularization in ill-posed linear problems. Malbet et al. [1] 
used this algorithm in a computationally expensive fashion since they were evaluating the 
gradient Vxg{Xn) numerically and also finding Vxg{Xn)^Vxg{Xn) at each iteration. The 
slower convergence rate due to the approximation in Eq. 11 is a minor drawback compared 
to the efficiencies gained by computing the approximate Hessian only. 

In practice, if it were possible to measure Ef perfectly with a precision better than 10~^, 
then only one series of numerical iteration would be sufficient to converge to a solution that 
would yield a contrast when the magnitude square of the errors in the estimate is below 
10~^^. However, in the presence of photon and camera noise, such a precision will never be 
achieved and a new estimate is needed each time the contrast improves. Moreover, since 
the estimation is based on DM diversity, and the actual DM shape is not fully known, the 
estimate of Ef becomes biased. As a consequence, while in theory the algorithm presented 
here could only consist of numerical iterations with only one measurement of Ef^ in reality 
it needs to include iterations with successive estimates of the field due to experimental 
limitations. 
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B. Electric Field Conjugation 

While some of the numerical difficulties of energy minimization could be alleviated by trying 
other quasi-Newton search techniques, such a path is not fruitfuU. Since we are minimizing 
the energy, we are not necessarily guaranteeing the contrast we are trying to achieve. That 
is, the line search approach never uses the fact that a coronagraph has been implemented 
to achieve high contrast. One alternative approach is the electric field conjugation (EFC) 
algorithm introduced by Give'on et al. [6]. In this algorithm, Give'on replaced the opti- 
mization with a root finding problem. Since we have a desired field in the image plane, 
Ed{^^v) — C{A{x^y)}^ the control problem becomes finding the DM settings required to 
achieve that field in <S, 



if we use a continuous version of the image plane, this problem is infinite dimensional and 
nonlinear in the a^. It is made tractable by discretizing the field at a finite number of points 
in the image plane. Give'on solved for the roots by assuming the DM surface height to be 
small and expanding the complex exponential in a Taylor series. 



where we have written the DM surface as the surface at the previous iterate, ipn^ pl^s a 
small correction, 

This lets us write the root finding problem at each iteration as a simple linear equation. 



Ef — E]j — 0. 



(16) 




(17) 



27rAo 



C{^/i(7A,X„)F}5X„+i = 



(18) 



where En = Ed - C{Ae"^''"'A''/^A+'A'''/'"} = Ed - {Ef)n is the error between the desired and 
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previously corrected aberrated electric field at the image plane and 5Xn-\-i is the change in 
the DM coefficients Xn at iteration n + 1. Since in all of the algorithms we are considering, 
we assume only an estimate of the electric field in <S, we can further simplify by expanding 
C{Ahnfk} and approximating it by it's zero order term, 

En + i'^C{AF}SXn+i = 0. (19) 

This is a linear equation in and can thus be solved for the DM increment given a 

previous measurement of the image plane field in S at step n, C{Ae^^^^^'~^'^''}. 

However, for it to be invertible, the same number of points in the image plane must be 
taken as there are coefficients, a^. Unfortunately, measurements are only available at the 
pixel spacings in the dark hole, which typically are much fewer in number than the available 
actuators. The result is an under-determined problem for the The approach taken 

in Give'on et al. [6] is to form the pseudo-inverse of C{AF}^ which is equivalent to finding 
the minimum norm solution for dX^^i. Unfortunately, this too can suffer from numerical 
difficulties. In particular, for a given coronagraph there is no guarantee that even the pseudo- 
inverse of C{AF} is well-behaved. This is particularly true for shaped pupils since many of 
the DM actuators are covered by the opaque regions of the mask. 
C. Stroke Minimization 

We solve many of these problems in our new approach, which we term "stroke minimization" . 
In stroke minimization, we eliminate the dimensionality problem by minimizing the finite 
number of coefficients in the DM surface expansion rather than the field itself. This has the 
effect of minimizing the average stroke of the DM actuators, an important criteria for the 
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small devices being used. We include in the minimization the constraint that the field meet 
the contrast requirement in the dark hole. 
The optimization problem is thus written, 

1 ^ 

minimize - ^ 

2 k=i 

subject to £s < 10"^ 

where £s — {Ef^ Ef)^ is the integrated intensity in the dark hole given by Eq. (6) and C is 
the contrast desired. 

While we have solved the dimensionality problem, this is still a difficult nonlinear pro- 
gram. A common solution approach is to expand the constraint in a Taylor series about the 
previous coefficient settings and replace the optimization with a quadratic subprogram at 
each step. Returning to our matrix notation where Xn is the matrix of coefficients at step n 
and 5Xn+i is the change in actuator settings, we can approximate the field at the DM, Eq^ 
inEq. (5), 



Eo^A{x,y) 



Xn 



(20) 



where the Jacobian of h is the row matrix, 



_ dh 



The integrated intensity can now be approximated using this first order expansion, 

Ss^ I I [C{AK} + C{AJr,}5Xr,+i\* [C{AK} + C{AJj(5X„+i] d^di^. (22) 



Multiplying through gives us, 



£s ^ J J^C{Ahn}*C{Ahn}dCdri+2ul^J J^C{Ahn}*C{AJ^}SXn+id^drj 



(£s)n 

14 



j jjXl^^C{AJr,yC{AJr,Y5Xr,+^d^dri 



where we have noted that the first term is just the measured energy in the dark hole at 
iteration n. This lets us write the quadratic subprogram, 

minimize ^(X^ + 5Xn^ifW-^{Xn + 
subject to 5Xl^-^MJXn^i + BJXn+i + dn < 10"^ 

where we have generalized to allow for an arbitrary weighting among the coefficients, through 
the matrix W, and, 



4 = J J^C{AKrC{AK}d^dr] 

^ I^C{AKrC{AJn}d^dr]^ 
Mn ^ j j^C{AJnyC{AJnYdidr]. 



Because 8s is always positive, this is a convex quadratic program and thus very efficient 
global solvers are available given estimates of d^, S^, and M^. One simple approach is to 
augment with Lagrange multiplier, /i, and solve the first order optimality condition as a 
function of /i. However, as with the other algorithms, we assume only estimates of the field 
in the image plane are available at each iteration, C{Ahn}' For the Jacobian, we proceed 
as presented previously by using a constant value obtained assuming ijjn = and ignoring 
the cross term between DM infiuence function and aberration. 



7 _ 

^'-dx 



= (23) 



This greatly reduces the computations needed in the algorithm since the state does not 
appear in the gradient. The augmented cost function can be written using the lagrange 
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multiplier /i, 

£m = l{Xr, + 5Xr,+ifW-\X^ + 5Xr,+i)+fi{SX^+,Mo5Xr,+i + BJXr,+, + d„-10-'') (24) 
The optimality condition corresponds to solving for the zeros of the derivative of Sm^ 

5Xr,+M = -{ill + 2WM^)-\-X^ + WBl) (25) 

/i 

where, 

5„ = 25R|y j^C{AKyC{AJo}did7i^ = {Ef,i^C{AF}^ (26) 
Mo = J I C{AJorC{AJo}^d^dri = -{^f {C{AF}, C{AF})s • (27) 

In order to find the optimal fi^^ we use a line search that finds the smallest fi — fi^ such 
that the contrast constraint is satisfied. Indeed, all the /i-dependent terms in the modified 
cost function correspond to a penalty that weighs the relative importance of contrast with 
respect to actuator minimization. Thus, the more stringent the contrast constraint, the 
larger fi^ becomes. Based on this qualitative approach, we start our algorithm with a small 
/io, compute X'*'(/io), simulate the propagation of these DM commands through the system, 
and increase /i until the contrast constraint is satisfied. We also begin with a smaller contrast 
target to ensure small DM changes at each iteration and slowly increase the target contrast 
until the goal of 10~^^ is reached. We call this algorithm "Stroke Minimization," since it finds 
the smallest deformations that achieve a target contrast. In Fig. 1 we show, via simulation, 
the PSF that results from one iteration of the Stroke Minimization algorithm at A = Aq, 
with Cxarget — 10~^^. Fig. 2 shows a comparison between the strokes obtained using a direct 
minimization of Ss and Stroke Minimization, for which the strokes are smaller by a factor 
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of two to five. For this simulation, the DM was modeled using an influence function basis, 
similar to the algorithms we use on our testbed. 
D. Experimental results 

1. High contrast PSFs and convergence curves 

The optical layout of the Princeton High Contrast Imaging Laboratory is shown in Fig. 3. 
It consists of two six-inch off-axis parabolic mirrors to coUimate and refocus the beam, a 
shaped pupil coronagraph, and two sequential DMs for control. The DM(s) are 1 cm on a 
side, and neither is located in a plane conjugate to the pupil. For the experiment presented 
here, we use only the image plane camera and control only one of the two DMs, while no 
voltage is applied to the other one. The dark hole obtained using the Stroke Minimization 
algorithm is shown in Fig. 4 and exhibits a dark hole 2 orders of magnitude deeper than the 
non-corrected one. The dark hour glass shape in the image plane corresponds to a binary 
mask which suppresses the bright central core and vertical wings of the PSF. Such an image 
plane mask is used to mitigate the limited dynamic range of our camera. 

2. Current Contrast Limitations at the Princeton Testbed 

Stroke Minimization also proves to be an effective diagnostic tool for the limitations of the 
testbed. Fig. 5 shows experimental convergence curves obtained using the Stroke Minimiza- 
tion algorithm. Each panel illustrates the contrast versus iteration curve for a given target 
contrast. We use the following notations: 

• ^Target' Target coutrast setting the contrast constraint of the wavefront control algo- 
rithm 
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• /^f': Integrated magnitude square of the estimated field in the Dark Hole at each 
iteration. 

• Integrated intensity actually measured in the Dark Hole at each iteration 

The last three panels of Fig. 5 illustrates this divergence. This divergence is due to the 
inherent bias of our DM diversity estimator. We note a discrepancy between the estimated 
intensity, /^^, and the actual intensity in the dark hole, I^^ . This implies an estimation 
bias that can be caused either by incoherent light landing in the dark hole or a systematic 
error in the estimation scheme. When Irarget — 1-6 x 10~^ the iterative loop diverges, 
primarily due to the large bias in the field estimate. Note that the contrast in the dark 
hole before correction in Fig. 5 is 10~^, an order of magnitude better than in Fig. 8. This 
estimate bias arises because the estimation algorithm relies upon a linearity assumption of 
the DM, perfect knowledge of the influence function, and a perfect knowledge of the voltage 
to deformation transfer function. When applied to the Princeton testbed, these assumptions 
produce an intrinsic bias in the estimated field. Each iteration of the stroke minimization 
algorithm minimizes the DM deformation under the constraint that the norm squared of the 
sum of estimated field and effect of the DM is below Irarget- If we decompose the estimate 
as 5Ef + Ef{Xn), where 5Ef is the bias, then the constraint can be written as: 

^Target (28) 
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where 5Ef is the estimation bias and || |p is the norm square. This imphes that the constraint 
in the Stroke Minimization loop is actually 

\\6Ef\\^ + \\Ef{Xr,) + C{AF}6Xr,\\^ - 2^ [{6Ef, Ef{Xr,) + J < Irarget (29) 

If we write the integrated value of the bias as Ii^ias = H^^^'/IP, the Cauchy Schwartz inequality 
yields: 

^Target -^bias + 2yfl^s\fSs (30) 

and the quadratic contrast constraint is more stringent than it is supposed to be. When 
hias < ^Target the algorithm is quite insensitive to the estimation error, as seen on the first 
three panels of Fig. 5. However, when these two quantities become similar then at some 
iteration rii the intensity constraint becomes too stringent; the algorithm seeks to correct for 
wavefront errors that are larger than the ones actually present in the testbed. This yields 
a set of deformation coefficients that is too large for the assumptions of the algorithm to 
be valid. The last three panels of Fig. 5 illustrates this divergence. The upwards trend 
in the last few iterations indicates that the algorithm might diverge if ran for a few more 
iterations. When tested the only the last panel actually diverges, for the two other case the 
algorithm lead to a oscillatory regime in contrast. Solutions to circumvent these limitations 
include estimation algorithms that do not use the DM as a source of diversity, or adaptive 
algorithms that build an on-the-fly model of the non-linear response of the DM and include 
that model in the estimation stage. 
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3. Symmetric Dark Hole with two DMs 

Another value of the stroke minimization algorithm is that it is easily modified to incorporate 
multiple DMs. We show in Pueyo & Kasdin [9] how multiple DMs are necessary to achieve 
symmetric dark holes on both sides of the image plane and for achieving high contrast 
in broader bands. The ability to provide a wavelength independent lever for wavefront 
correction is the main advantage of controllers based on two sequential DMs since it will 
ultimately enable broadband observations of exo-planets and thus greatly facilitate their 
spectral characterization. The actual implementation of control algorithms in broadband 
using these methods can be either done using a series of monochromatic estimations, as 
shown by Give'on et al. [6], or a single wavelength estimation coupled with some priors 
on the symmetries of the PSF. In either case, once the estimation is complete the stroke 
minimization algorithm can be applied in order to retrieve the DM commands. Discussing 
relative performances of wavefront retrieval methods under polychromatic light is beyond the 
scope of this paper, and here we chose to only focus on the intricacies of the implementation 
of a two DMs control algorithm. As a consequence we will only present on monochromatic 
results, using the same estimator as in § 2D, using the two DMs to produce a symmetric 
dark hole. 

A. Monochromatic Stroke Minimization with two sequential DMs: general algorithm 
Here we do not delve into to the details of the several single DM correction algorithms 
presented above and solely focus on Stroke Minimization. To begin, we write the matrix 
of actuator commands as a concatenation of the coefficients of each DM, X = [X^-^^X^^^]. 
When we consider the case where DM2 is conjugated with the shaped pupil (or the Lyot 
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plane of the coronagraph) and separated from DM1 by a distance the field at DM2, before 
the final Fourier Transform, can be written as 

EDM2ix,y) = Aix,y)h^^\jx,X) (31) 

where 

/i(^)(7a,X) = e^^(^'^)+^^^^(^'^)j;[e^^'^^'^(^'^)]e^^'^^^^ (32) 

'0^^) and V^^^) stand respectively for the surface of DM1 and DM2 at the nth iteration, is 
the Fresnel propagation between two surfaces separated by a distance z, and A{x^y) is the 
pupil apodisation. Using the coronagraph operator C notation, the field in the final image 
plane is then: 

The general form for the intensity resulting from the effects of the two DMs is then 

V V 



J J^AX^C{AJn}*C{AJnrSXd^dr]. 



where, by virtue of the Unearity of the operator C, 

Jn^^\x„. (34) 

Thus, just as we presented above for the case of a single DM, we are seeking to solve the 
following optimization problem: 

1 ^ 

minimize - ^ a\ 

2 k=i 

subject to Es < 10"^ 
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We proceed using the same sub-quadratic programming approach and solve at each iteration 
the following subprogram: 

minimize ^(X^ + 5Xn+ifW-^{Xn + 
subject to 5Xl^-^MJXn+i + BJXn+i + dn < 10"^ 

However for the case of multiple DMs we need to reduce the modeling of the wavefront 
controller in such a way that M^, and dn can be directly computed either from wave- 
front estimate or the design parameters of the optical set up. Once again this is done by 
approximating by a constant value Jq^ which slows down the convergence rate of the al- 
gorithm but circumvents the high computational cost associated with the evaluation of the 
sensitivity matrix at each iteration. Next we present the model reduction we implemented 
for our experimental validation. 

B. Monochromatic Stroke Minimization with two sequential DMs: model reduction 
While dividing the problem of finding optimal DM strokes in a a series of quadratic sub- 
programs is a general method applicable to all types of wavefront control architectures, here 
we are interested in reducing the sub-program in such a way that: 

• The sensitivity matrix is computed only once, before any correction 

• The sensitivity matrix provides two degrees of freedom to correct on both sides of the 
image plane 

• The sensitivity matrix provides two degrees of freedom for wavelength independent 
wavefront errors and those proportional to 1/ A [9] 
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We show here how to reduce Eq. 32 in such a fashion. The main difference between Eq. 5 
and Eq. 32 is the presence of a Fresnel propagation between the two DMs . This propagation 
is what provides both the symmetric and wavelength levers. We start by calculating the 
jacobian of the effect of the DM in the plane where the aberrations are estimated, 

= i^e"x{x,y)+i'-^Mx,y) (35) 
A 

where F^^\x^y) is the matrix representing the basis function for the j th DM. We proceed 
as previously, namely we choose to evaluate this jacobian around X — Q and ignore the 
cross talk between the DMs and the aberrations. We also ignore the cross talk between the 
two DMs. These approximations lead to a Jacobian that is not exact but making them only 
lowers the convergence rate and does not change the final solution. The contribution of the 
second DM becomes: 

A A 

(36) 

In order to calculate the contribution of the first DM, that is not conjugated with the plane 
of the aberrations, we assume, as it is the case for shaped pupils, that the operator C is a 
fourier transform. Moreover we work under the angular spectrum approximation, and thus 
the impact of an out of pupil optics is only to multiply the electrical field in the image plane 
by a quadratic phase factor. The contribution of the out-of-pupil plane DM is then: 

C{iA{x, ^)^e"^(^'^)+^x^A(-.^)e^x^i>'^(-.^)jr^[e^¥^^'^(-.2/)]F(i)(a;, y)} (37) 
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where again we have used the fact that the gradient is computed around X = with no 
aberration. Thus the reduced and hnearized effect of the two DMs becomes 



C{AJn} ~ C{AJo} = i 



[C{a4'^} C{AJm 



(2)1 



(38) 



In our laboratory we have implemented the stroke minimization algorithm with dual DMs 
using this reduced Jacobian. Namely, at each iteration we solve the quadratic subprogram 

minimize ^{X^ + 5Xn+ifW-\Xn + 
subject to 5X^^^M5Xn+i + BJXn+i + 4 < 10"^ 



where 



dn = (C{Ah(^^'C{Ah(^^ 



C{Ah(^n, C{A4'^})^ {C{Ah(^n, C{AJr )}/^ 



M = 



'{a4\ c{a4'^})^ {c{a4% c{a4'^] 

'{a4% C{a4'^})^ {C{a4\ C{Ajf ] 
C. Stroke Minimization with two sequential DMs: symmetries and broadband lever 



Note that for small spatial frequencies the angular spectrum is small, ^(^^ + r/^) <C 1, 
resulting in the following first order Taylor expansion 



C{AJo} 



^^{e + 77')C{AFW} + i^C{AF('} i^C{AF(2)} 



(39) 



If we rearrange the DM commands as X = [Xi X2] = [Xi Xi + X2] then the sensitivity 
matrix becomes: 



C{AJo} ^ 



^^ie + v'')C{AFW} i^(^c{AF^^} + C{AF(2)}) 



(40) 
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Because F^^\x^ y) and A{x^ y) are real function in the pupil plane, C{AF^^^} and C{AF^'^^} 
feature a Hermitian symmetry in the image plane when C{ } is a fourier transform. Thus 
the first block of C{AJq} provides a lever that corrects for Hermitians field distributions 
in the image plane while the second one provides a lever that corrects for anti-Hermitian 
distributions. It is the combination of these two independent Hermitian and anti-Hermitian 
levers that allows us to create symmetric monochromatic dark holes. Note that in the case 
of a single DM correction there is only an anti-Hermitian lever, thus both sides of the image 
plane cannot be corrected independently, which constraints the controllable area to only 
one side of the PSF. These considerations are developed in Appendix A using a phasor 
representation. 

Moreover, each of these two independent levers has a different chromatic behavior: the 
hermitian one scales as and the anti-hermitian one scales at 1/A. In terms of wavefront 
errors this means that a two sequential DM controller can correct over a broadband for 
what are commonly called amplitude and phase errors in the pupil plane. However, im- 
plementing a broadband stroke minimization control algorithm requires a more elaborate 
estimate of the field in the image plane. One option is to obtain an estimate across the 
bandpass that provides a set of and dn for each wavelength so that the constraint of 
the quadratic subprogram becomes the integrated broadband intensity. An alternative is to 
force a monochromatic quadratic subprogram to use the Hermitian lever only to correct for 
the component of the estimated wavefront and to use the anti-Hermitian lever for its 1/A 
component. While critical for the feasibility of broadband wavefront control and thus the 
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detection of exo-planets, the development and implementation of such estimation algorithms 
is beyond the scope of this paper and will be presented in a future communication. For the 
remainder of this article we focus on experimental results that feature a monochromatic 
symmetric dark hole. 
D. Experimental results 

Fig. 6 shows the DM surfaces obtained using this algorithm to create the symmetric 10~^^ 
monochomatic dark hole that is presented on the top panel of Fig. 7. For these numerical 
simulations we have used square DMs of size D — ?> cm that are separated by z = 1 m. In 
this section, we present the first results of a symmetric dark hole using two DMs in sequence 
to correct for errors on both sides of the image plane. This experiment was performed in 
monochromatic 635 nm light using the stroke minimization algorithm as described in § 3 A. 
As before, the estimate of the wavefront was obtained using an algorithm based on the 
application of diversity on the surface of one of the DMs. As shown in Fig. 3, neither of 
the two DMs is in a plane conjugate to the shaped pupil. The propagation from each DM 
to the pupil plane is taken into account using the angular spectrum approximation, just as 
shown for the first term of Eq. 39 

Figure 8 shows the aberrated image prior to correction as well as the image after 60 
iterations of the stroke minimization correction algorithm. In addition, the figure shows a 
contrast plot as a function of iteration. The Dark Hole is from 7-10 \/D mx and -3 to 3 A/D 
in y. The average contrast between the two sides of the image plane before any correction 
is at 1.2 X 10~^ with the right side starting out worse than the left side. After 60 iterations, 
the contrast on both sides of the image has reached to 2.5 x 10~^. This figure shows that the 
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stroke minimization algorithm allowed us to improve the on-axis light extinction by almost 
two orders of magnitude. 

4. Conclusion 

In this paper we presented a novel general method to solve the non-linear inversion problem 
associated with the correction of quasi-static wavefront errors. This novel algorithm that we 
named Stroke Minimization circumvents the dimensionality of the problem and allows a se- 
lection of the regularization parameters that is directly related to the target contrast desired 
in the Dark Hole, where exo-planets are expected to be seen. It can also easily be generalized 
to multi-DM systems. In this communication, we used this algorithm to accomplish the first 
experimental proof of a symmetric high contrast PSF obtained using two sequential DMs. 
This is a significant experimental milestone for the field of high contrast imaging since it not 
only doubles the search space of coronagraphs, but also proves that amplitude and phase 
errors can be simultaneously corrected. A full chromatic characterization of this solution 
will be presented in a subsequent communication. 

Appendix A: Physics of wavefront correction 

We present here a qualitative explanation of the physics involved in a single-DM wavefront 
compensator. We earlier mentioned that such a device could only create a dark hole in 
half of the image plane (c.f.. Brown & Burrows [2]), and we show here how it achieves such 
a feature. First, consider a phase error in the pupil plane that is composed of only one 
harmonic component: 

y) = e'^ cos(^(m«)+0) ^^^^ 
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Then, it can be corrected using a deformable mirror to conjugate the field in the pupil plane: 

E^Zi^. y) = e-'^^"<^("^"+"^)+^) (A2) 

This correction is perfect for all the wavelengths. Now consider the case of an amplitude 
error: 

Eabbi^^ y) = COS (^^{mx + ny) + (A3) 
and of a DM surface that is such that, to first order: 

^^M(2;,y) = -iysm(^(m2; + ny) + 0) (A4) 

Then the residual field in the pupil plane is: 

^fleT(^> y) = 1(1 - A)e<^(--+"2/)+^) + ^(1 + A)e-<^(--+"^)+^) (A5) 
2 Aq 2 Aq 

When A = Aq then the positive spatial frequencies are corrected while the negative ones are 
not compensated. We illustrate this feature using a phasor representation. Fig. 9 shows how 
the spatial variations of amplitude and phase errors are represented as pulsating phasors that 
can be decomposed into two rotating phasors in the complex plane. The clockwise rotating 
phasor corresponds to the contribution of the ripple in the right half of the image plane and 
the anti-clockwise in the left half of the image plane. One can cancel the clockwise component 
of an amplitude error using the anti-clockwise component of a phase deformation introduced 
by a DM. This concept is illustrated in Fig. 9. This is exactly the approach carried out in the 
"speckle- nulling" algorithm (ref Borde Traub), where there is no wavefront estimation and 
the alignment of the phasor occurs via a trial and error process which considerably lengthens 
the convergence time (See § 1). A similar analysis can be performed for wavefront correction 
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using 2 DMs showing that the wavefront actuator can correct both amphtude and phase 
aberrations under a broadband illumination. This capability was first shown in Shaklan 
et al. [8]. Using the set up of Fig. 3, where DM2 is in a plane conjugate to the final imaging 
lens. DM2 can be used to correct phase errors as shown by Eq. A2. Amplitude errors 
can then be compensated using a combination of DM1 and DM2. Assume the amplitude 
aberration, located at DM2, is such that: 

Epup,abb{x, y) = COS {^{mx + ny) + (A6) 

We choose the surface of DM1, such that the linear contribution of DM1 to the field is: 

Edmi,pup{x, y) = — cos f ^(mx + ny) + (A7) 

Then, using the results of [9] and [10], the propagation from DM1 to DM2 of this field is 
given by: 

Edmi,pup{oo, y) = -i- . . ^ , ' cos -—{mx + ny) + (A8) 

We work in the low-to-mid spatial frequency regime, so here again we can assume that 
^^^^^2^^ ^ ^ 1- With this first order approximation of the angular spectrum factor the 
contribution of DM1 at DM2 becomes: 

, , .Ao /27r, , , , [271 ^ 

Et 



^WMi,pup{x, y) = ^^;^^(-^2 _^^2-) (^(^^ + ^y) + ^) - COS (^-^{mx + ny) + 

(A9) 

Therefore choosing: 

Edm2,pup{x, y) = cos (^(mx + ny) + ^) (AlO) 
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yields a broadband cancellation of the amplitude error. These considerations are illustrated 
using a complex phasor representation on Fig. 10, where DM1 corrects for the amplitude 
errors and DM2 for the phase errors. Note that there is an important assumption underlying 
this result: the small angular spectrum regime is required in order to obtain the broadband 
property of the two DMs controller. This yields an outer working angle limit, dependent on 
the optical design of the controller, that has been derived in Pueyo & Kasdin [9]. 
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Fig. 1. Numerical results of the Stroke Minimization algorith using a 10 Shaped Pupil. 
Top: DM deformation in radians. Bottom: Log(Corrected PSF). 
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strokes Minimization 




log(Extinction) = - log(Contrast) 



Fig. 2. Comparison between the peak to valley actuator strokes necessary to correct a given 
wavefront error using Energy Minimization and Stroke Minimization. 
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Fig. 3. Optical layout of the princeton high contrast imaging testbed. For the experiment 
presented in § 2D, only one DM is used for wavefront control and only the image plane 
camera is used for wavefront sensing. 
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Aberrated Image Contrast Corrected Image Contrast 




Sky angle irA(3/D Sky angle IrAc/D 

Fig. 4. Aberrated (Left) and corrected (Right) PSF on the Princeton testbed in 
log (contrast). The wavefront is flattened so that half of the image plane exhibits a dark hole 
over a specifled region, in this case X = 7-10 \/ D and Y = -2.5-2.5 \/ D. This monochro- 
matic experiment used an illumination wavelength of 635 nm 
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Fig. 5. Experimental results for six different target contrasts. The top two curves correspond 

respectively to the maximum of the intensity in the dark hole and the maximum of the 

estimated intensity in the dark hole. The bottom two correspond to the average intensity 

in the dark hole and the average estimated intensity in the dark hole. Note that when the 

algorithm converges the average estimated intensity is equal to the target contrast. Also, 

note that these results were obtained on the Princeton testbed prior to the installation of 

the second DM, and therefore, the contrast limit is slightly better than other results shown 
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(a) Wavefront deform ation in radians ^,-1 (b) Wavefront deformation in radians 




Fig. 6. DM surfaces in radians obtained using the two DMs Stroke Minimization algorithm. 
The algorithm used here is designed to operate monochromatically and does not take ad- 
vantage of the broadband capabilities of the wavefront controller. Pupil size: D = 3 cm and 
DM separation: z = 1 m 
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Ao/D 

Fig. 7. Monochromatic PSF resulting from two DM wavefront correction using a monchro- 
matic Stroke Minimization algorithm. D = 3 cm and z = 1 m. 
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(a) Aberrated Image Contrast 
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(b) Corrected Image Contrast 
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(C) Contrast Plot 
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Fig. 8. Aberrated image (left) and corrected image(right) of the 2 DM stroke minimization 
symmetric dark hole experiment. Bottom: A plot of contrast vs. iteration in each of the 



two dark holes and in the combination of the two. 
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Fig. 9. Half dark hole correction using one DM. Left column: a phase aberration can 
theoretically be compensated at all wavelengths by a matching DM setting to cancel out the 
total electric field on both sides of the optical axis. Right column: In the case of amplitude 
aberrations, a DM setting can be found that will exactly compensate the amplitude error 
on one side of the optical axis, and at a single wavelength 
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Amplitude correction with DM1 Phase correction with DM2 




Fig. 10. Broadband amplitude correction using 2 DMs: complex phasor illustration. Left 

Column: one can find a DMlsetting that will cancel the amplitude error achromatically on 

both sides of the optical axis after propagation through the system (equivalent to a phasor 

rotation in the angular spectrum approximation). Right column: the phase error induced 

by DM1, together with the phase error accumulated through the system, is then taken out 

at all wavelengths by thesecond DM, resulting in broad-band 2-sided light cancellation. 

41 



